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We study the influence of astrophysical formation scenarios on the precessional dynamics of spin- 
ning black-hole binaries by the time they enter the observational window of second- and third- 
generation gravitational- wave detectors, such as Advanced LIGO/Virgo, LIGO-India, KAGRA and 
the Einstein Telescope. Under the plausible assumption that tidal interactions are efficient at align- 
ing the spins of few-solar mass black-hole progenitors with the orbital angular momentum, we find 
that black-hole spins should be expected to preferentially lie in a plane when they become de- 
tectable by gravitational-wave interferometers. This "resonant plane" is identified by the conditions 
A$ = 0° or A$ = ±180°, where A$ is the angle between the components of the black-hole spins 
in the plane orthogonal to the orbital angular momentum. If the angles A$ can be accurately mea- 
sured for a large sample of gravitational-wave detections, their distribution will constrain models 
of compact binary formation. In particular, it will tell us whether tidal interactions are efficient 
and whether a mechanism such as mass transfer, stellar winds, or supernovae can induce a mass- 
ratio reversal (so that the heavier black hole is produced by the initially lighter stellar progenitor). 
Therefore our model offers a concrete observational link between gravitational-wave measurements 
and astrophysics. We also hope that it will stimulate further studies of precessional dynamics, 
gravitational-wave template placement and parameter estimation for binaries locked in the resonant 
plane. 

PACS numbers: 04.25.dg, 04.70.Bw, 04.30.-w 



I. INTRODUCTION 

The inspiral and merger of stellar-mass black-hole 
(BH) binaries is one of the main targets of the future 
network of second-generation gravitational-wave (GW) 
interferometers (including Advanced LIGO/Virgo [T], 
LIGO-India [2] and KAGRA fS^) and of third-generation 
interferometers, such as the proposed Einstein Telescope 
d] . Typical GW signals from these binaries are expected 
to have low signal-to-noise ratios, and must therefore be 
extracted by matched filtering, which consists of comput- 
ing the cross-correlation between the noisy detector out- 
put and a predicted theoretical waveform, or template 
(see e.g. [5]). The number of observationally distinguish- 
able merger signals should be extremely large, both be- 
cause of the large and strongly mass-dependent number 
of cycles in each signal and because the emitted waveform 
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depends sensitively on as many as 17 different parame- 
ters, in the general case where the BHs are spinning and 
in eccentric orbits. The difficult task of exploring such 
a high-dimensional space can be simplified if nature pro- 
vides physical mechanisms that cause astrophysical bi- 
naries to cluster in restricted portions of the parameter 
space. 

In this paper we consider one mechanism to preferen- 
tially populate certain regions of parameter space: the 
post-Newtonian (PN) spin-orbit resonances first discov- 
ered by Schnittman [6^. Unfortunately, very few of the 
existing population-synthesis models of compact-binary 
formation (see e.g. [3 [5]) include self-consistent predic- 
tions for BH spins. To highlight the significance of spin- 
orbit misalignment and resonances, we adopt a simpli- 
fied model for binary BH formation. We use this model 
to generate initial conditions for our compact binaries, 
and then integrate the PN equations of motion forward 
in time using an extension of the code used by some of 
us in previous studies of supcrmassive BH binaries [51- 
[TT] . Our analytically tractable model captures (at least 
qualitatively) many of the detailed physical effects influ- 
encing the evolution of BH spins. Within this framework 
we carry out Monte Carlo simulations to study the statis- 
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Merging BH binary 




Resonant plane locking 
sinA$ — > (equilibrium) 



Free precession 
sinA$ -> ±1 (pile-up) 




Spin alignment 
A* 0° 

0,2 ^ 0° 



Spin anti-alignment 

A$ ±180° 
012 -^01 + 02 (tail) 



FIG. 1. Schematic summary of our predictions for the spin orientation of BH binaries as they enter the LIGO/Virgo band. If 
tides efficiently align the spin of the secondary with the orbital angular momentum prior to the second supernova, resonant-plane 
locking will drive sinA$ — >■ 0, while in the absence of tides the spins will precess freely, piling up around sinA<I> — > ±1 near 
merger. When tides are efficient, if the primary star evolves into the less massive BH (reversed mass ratio) the PN evolution 
will drive A$ — ^ 0°,6i2 — > 0°. If instead the primary star evolves into the more massive BH (standard mass ratio) the PN 
evolution will drive A$ — > ±180°, 9i2 -^01+02, generating a tail in the distribution of 612 to larger values. See Eqs. ([2| and 
S for definitions of these angles. 



tical distribution of BH spins when they enter the GW- 
detection band of second- and third-generation detectors. 

Before summarizing our results, we first introduce 
some notation. Consider a BH binary with component 
masses mi > TO2, total mass M = rrii + m2 and mass 
ratio q = 'm2/mi < 1. The spin Si of each BH can be 
written as 



Si = Xi 



-Si, 



(1) 



where < < 1 (i = 1, 2) is the dimensionless spin mag- 
nitude and a hat denotes a unit vector. Our goal is not 
to rival the complexity of existing population-synthesis 
models of compact-binary formation, but rather to inves- 
tigate specifically those astrophysical ingredients which 
affect the spin dynamics. We therefore focus on maxi- 
mally spinning BH binaries with mass ratio q = 0.8, a 
typical value predicted by population-synthesis studies 
(cf. e.g. Fig. 9 of mi). 

Let us define 6i to be the angle between each spin Si 
and the orbital angular momentum of the binary L, 9i2 
to be the angle between Si and S2, and A$ to be the 
angle between the projection of the spins on the orbital 
plane: 



cos6'i = Si • L, cos 02 = S2 • L, 



(2) 



COS01, = 81-82, cosA<i>=-El-^.-E^-^. (3) 

|Si X L| IS2 X L| 

As we will demonstrate below, the physical mechanisms 
leading to the formation of the BH binary leave a char- 
acteristic imprint on the angles A$ and ^12. This has 
implications for GW data analysis and, even more strik- 
ingly, for GW astronomy: at least in principle, measure- 
ments of spin orientation with future GW detections can 
constrain the astrophysical evolutionary processes that 
lead the binary to merger. 

All BH binaries with misaligned spins {6i 7^ 0) experi- 
ence PN spin precession as they inspiral towards merger. 



Although ensembles of BH binaries with isotropic spin 
distributions retain their isotropic distributions as they 
inspiral |T3j, anisotropic spin distributions can be sub- 
stantially affected by PN spin precession [5]. In particu- 
lar, binaries can be attracted towards PN spin-orbit res- 
onances in which the BH spins and orbital angular mo- 
mentum jointly precess in a common plane ("resonant- 
plane locking" ) . Binaries in which the two BH spins and 
the orbital angular momentum do not share a common 
plane at the end of the inspiral are said to precess freely. 
Binaries can become locked into resonance if they satisfy 
the following conditions at large separations: 

i) comparable but not equal masses (0.4 < q 7^ 1), 

ii) sufficiently large spin magnitudes {xi ^ 0.5), 

iii) unequal spin misalignments {9i 7^ 02)- 

If these conditions are satisfied, the spin distribution 
of an ensemble of binaries will be strongly influenced by 
the PN resonances although every individual member of 
the ensemble will not necessarily become locked into res- 
onance. In ensembles of binaries for which 9i < 62 at 
large separations, the two spins tend to align with each 
other, so that A$ -> 0°, 6*12 ^0°. If instead Oi > 62, 
the projections of the BH spins on the orbital plane tend 
to anti-align, so that A$ 180°, 6*12 9i + 6*2. The 
mass ratios for which resonant-plane locking is effective, 
given by condition i) above, are typical for the stellar- 
mass BH binaries detectable by Advanced LIGO/Virgo 
(cf. Fig. 9 of [12). The spin magnitudes Xi of newly 
formed BHs are highly uncertain, but observations of ac- 
creting BHs in binary systems indicate that their spins 
span the whole range < < 1 allowed by general rela- 
tivity [Mj [15] . Many BH-BH systems may therefore sat- 
isfy condition ii) above. In contrast, we would not expect 
resonance locking in binaries in which one or both mem- 
bers are neutron stars, as they are expected to have small 
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spinf[^ Whether the spin misahgnments of the ensemble 
of BH binaries detectable by Advanced LIGO/Virgo is 
asymmetric (satisfying condition iii) above) is a primary 
consideration of this paper. 

Astrophysical formation channels determine the initial 
conditions for PN evolutions in the late inspiral. As a re- 
sult they determine whether resonant locking can occur, 
and which resonant configuration is favored. Here we 
introduce a model for BH binary formation that allows 
us to establish a link between binary-formation channels 
and the near-merger spin configurations of processing BH 
binaries. 



A. Executive summary 

Our main findings are summarized schematically in 
Fig. [l] Supernova (SN) kicks tilt the orbit, producing 
a misalignment between the orbital angular momentum 
and the orientation of the spins of the binary members 
[25| . As a result, the main factors determining the spin 
alignment of a BH binary are the magnitude of SN kicks 
and the possibility that other physical effects may realign 
the spins with the orbital angular momentum in between 
SN events. Dominant among these physical effects (aside 
from the SN kick itself) are the efficiency of tidal inter- 
actions and the possibility of a mass-ratio reversal due 
to mass transfer from the initially more massive, faster 
evolving progenitor. 

Tides affect the binary in two significant ways: they 
align the spins of stellar BH progenitors with the or- 
bital angular momentum and they reduce the binary ec- 
centricity. Additionally, tides force stars to rotate syn- 
chronously with the orbit, increasing the likelihood of a 
large BH spin at collapse and implying that our results 
will depend only mildly (if at all) on the initial stellar 
spin. Consider the evolution of the system between the 
two SN events, when the binary consists of a BH and 
a non-degenerate star. If tidal interactions are efficient 



(a reasonable assumption, as we argue in Appendix A 6 ) 
they tend to align the star (but not the BH) with the or- 
bital angular momentum. This introduces an asymmetry 
in the angles {6i, 62) which is critical to determining the 
spin configuration at the end of the inspiral. 

Mass transfer can change the mass ratio of interact- 
ing binaries. Since the main-sequence lifetime of a star 



^ Relativistic calculations of neutron star structure suggest that 
Xi S 0-7 for uniform rotation and physically motivated equations 
of state |16H18I , but the spin magnitudes of neutron stars in 
binaries observable by Advanced LIGO are likely to be much 
smaller than this theoretical upper bound 1191 I20| . The spin 
period of isolated neutron stars at birth should be in the range 
10-140 ms 1211 . or Xi 0.04. Accretion from a binary companion 
can spin up neutron stars but is unlikely to produce periods less 
than 1 ms, i.e. Xi $ 0-4 1221 . The fastest-spinning observed 
pulsar has a period of 1.4 ms, {xi ~ 0.3) 1231 : the fastest knowrn 
pulsar in a neutron star-neutron star system, J0737-3039A, has 
a period of 22.70 ms (xi ~ 0.02) [2l]. 



is a decreasing function of its mass, the initially more 
massive star in a binary is expected to collapse first. If 
mass transfer from this star to its less massive compan- 
ion is insufficient, which we will refer to as the standard 
mass ratio (SMR) scenario, the initially more massive 
star will go on to form the more massive member of the 
BH binary. We cannot however rule out the possibility 
that prior to the first SN, the initially more massive star 
overfiows its Roche lobe and donates mass to its initially 
lighter, longer-lived companion. This mass transfer may 
produce a mass-ratio reversal, so that the heavier BH in 
the binary forms second: we will call this the reversed 
mass ratio (RMR) scenario. According to population- 
synthesis models, mass-ratio reversal happens for a siz- 
able fraction (typically from ^ 10% to 50%) of the total 



number of BH binaries (cf. [12J and Table III below) 



Since BHs are relatively immune to the effects of tides, 
the spin of the first BH to form will be more misaligned 
than the spin of the second BH, as this misalignment will 
have accumulated due to the kicks generated during both 
SN events. Therefore, in the SMR scenario BH binaries 
will have 9i > 02 at formation, and thus A$ ~ ±180° 
by the time they enter the GW-detection band. On the 
other hand, in the RMR scenario BH binaries initially 
have 9i < 82, so that by late in the inspiral A$ ~ ±0°, 
and furthermore the spins are nearly aligned with each 
other (i.e., 9i2 ~ 0). In summary, whenever tidal inter- 
actions are efficient, our model predicts that BH spins 
should preferentially lie in a "resonant plane" (identified 
by the conditions A<I> — 0° in the RMR scenario, and 
A$ — ±180° in the SMR scenario) when they become 
detectable by GW interferometers. 

A third (more unlikely) possibility is that tidal inter- 
actions are not efRcient. In this case, binaries form with 
di ~ 02 and will not become locked into resonant configu- 
rations. Our simulations show that binaries will preferen- 
tially have A<I> ~ ±90° . Because the most likely values of 
A$ in the three scenarios (RMR, SMR and no tides) are 
mutually exclusive, GW measurements of a statistically 
significant sample of values of A$ will provide important 
astrophysical information on compact-binary formation 
scenarios. In particular, they will tell us whether tidal 
interactions are efficient, and (if so) whether mass trans- 
fer can produce mass-ratio reversals. 

Fig. [2] makes these conclusions more quantitative by 
showing three histograms of A$ (left) and 6*12 (right), 
corresponding to snapshots taken at different times dur- 
ing the inspiral. The distribution of A<1> is flat at large 
separations (dotted lines, corresponding to early times 
and small orbital frequency) because spin-spin couplings 
are weak, and the BH spins simply precess about the 
orbital angular momentum. If tidal alignment is effi- 
cient, in the late inspiral the BH spins lock into equi- 
librium configurations with either A<I> = 0° or A<I> — 
±180°. This effect is clearly visible at GW frequencies 
/gw = 1 Hz, roughly corresponding to the lowest cutoff 
frequency of third-generation detectors like ET, and it is 
even more pronounced when the binaries enter the Ad- 



4 




FIG. 2. (Color online.) Left: Probability distribution of the angle between the projections of the spins on the orbital plane 
A$. As the binaries inspiral, the GW frequency /gw increases from 0.01 Hz (dotted blue lines) to 1 Hz (dashed red lines) and 
later 20 Hz (solid black lines). Under the effect of tides the PN evolution brings the spins in the same plane (A<E> — >■ 0°, ±180°), 
both in a reversed mass ratio (RMR, top panel) and in a standard mass ratio (SMR, middle panel) scenario. When tidal effects 
are removed (bottom panel, where we show both RMR and SMR binaries) the spins precess freely and pile up at A$ — ±90°. 
Right: Probability distribution of the angle between the two spins 9i2. In the RMR scenario (top panel) the spins end up 
almost completely aligned with each other, i.e. most binaries have O12 — 0°. In the SMR scenario (middle panel) and in the 
absence of tides (bottom panel, where again we show both RMR and SMR binaries) a long tail at large values of 6^12 remains 
even in the late inspiral. All simulations shown in this figure assume that kick directions are isotropically distributed. Error 
bars are computed assuming statistical Poisson noise. 



vanced LIGO/Virgo band at /gw — 20 Hz. If tides are 
artificially removed, free precession during the late stages 
of the inspiral slows down the evolution of A$ when the 
components of the spin orthogonal to the orbital angu- 
lar momentum are also orthogonal to each other, causing 
binaries that are not locked into resonance to pile up at 
A$ = ±90°. 

Let us stress again that the statistical effect of res- 
onances is clearly visible at /gw = 20 Hz, i.e. when 
BH binaries enter the Advanced LIGO/Virgo band. GW 
measurements of A$ can therefore be used to constrain 
uncertainties in BH binary-formation scenarios. The in- 
clusion of resonant effects in population-synthesis models 
(combined with a statistically significant sample of GW 
measurements of A<I>) has the potential to constrain var- 
ious aspects of the models, such as the efficiency of tides, 
stable mass transfer, common-envelope (CE) evolution, 
SN kick velocities, and the metallicity of BH progenitors. 



B. Outline of the paper 

The rest of the paper provides details of our as- 
trophysical model and a more detailed discussion of 



the results. In Section |TT] we introduce our fiducial 
BH binary-formation channels, which are based on de- 
tailed population-synthesis models, as described at much 
greater length in AppendixjX] In order to focus on spin ef- 
fects, we fix the component masses to two representative 
values. We assume that SN kicks follow a Maxwellian 
distribution in magnitude. We also assume that the 
kicks are distributed in a double cone of opening angle Ob 
about the spin of the exploding star and, to bracket un- 
certainties, we consider two extreme scenarios: isotropic 
(6*6 = 90°) or polar (6*6 = 10°) kicks. 

Section |III| summarizes the results of evolving these 
BH binaries under the effect of gravitational radiation 
down to a final separation of lOGM/c^. We demonstrate 
that spin-orbit resonances have a significant impact on 
the observable properties of our fiducial BH binaries. Al- 
though we have only explored a handful of evolutionary 
channels and component masses, in Section [rv| we argue 
that the scenarios described in Fig. [T] are broadly appli- 
cable: kicks, tides, and the mass-ratio distribution con- 
trol spin alignment. We explore the sensitivity of these 
three features (and hence of the observable distribution of 
resonantly-locked binaries) to several poorly constrained 
physical inputs to binary-evolution models, and we argue 
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that GW observations of precession angles could provide 
significant constraints on binary formation channels. Fi- 
nally, in Section |V] we describe the implications of our re- 
sults for future efforts in binary-evolution modeling and 
GW detection. 

To complement and justify the simple astrophysical 
model proposed in Section |lTj in Appendix [K\ we describe 
in detail the rationale underlying the model and its re- 
lationship to our current understanding of binary evo- 
lution. Appendix [X] should provide a useful resource to 
implement (and possibly improve) the Monte-Carlo algo- 
rithm described in the main text. 



II. ASTROPHYSICAL MODEL OF THE 
INITIAL CONDITIONS FOR SPIN EVOLUTION 

Isolated BH binaries do not emit electromagnetically 
and hence have yet to be observed. Despite this lack 
of evidence, they are a likely outcome of the evolution 
of massive stellar binaries. The rate at which they form 
can be inferred from observations of their progenitors and 
systems like binary neutron stars that have similar for- 
mation channels. Formation rates can also be calculated 
theoretically using population-synthesis models such as 
Star Track [^BH25] , which builds upon previous ana- 
lytical studies of single [29j and binary stellar evolution 
[30|. 

Most studies of compact-binary formation do not keep 
track of the magnitude and orientation of BH spins, 
and those that do (see e.g. O |H1 [31]) neglect general- 
relativistic effects in the late-time evolution of the bi- 
nary. One of the goals of our study is to fill this gap. For 
example, the version of the StarTrack code used in [7] 
assumed that both Si and S2 remained aligned with the 
initial direction of the orbital angular momentum L. The 
evolution of L itself was performed by applying energy 
and angular-momentum conservation when compact ob- 
jects are formed (and kicked) as a result of gravitational 
collapse. This approach is suitable for binaries in non- 
relativistic orbits, like observed X-ray binaries [8l l3T] . 
but it may not be appropriate for merging binaries, that 
are interesting both as GW sources and as progenitors of 
short gamma-ray bursts [7]. Since existing BH binary- 
formation models preserve the mutual alignment of the 
spins with the initial direction of L, all BH-BH binaries 
are formed with 61 = 62- Later models of mixed BH 
X-ray binaries do allow for the possibility of asymmetric 
spin configurations via accretion f8], but to the best of 
our knowledge no such studies have been published for 
the BH-BH case. Since PN resonance locking only occurs 
when 9i 02, its effects are excluded by construction in 
the BH binary models available in the literature. 

Here we develop a slightly more complex (and pre- 
sumably more realistic) model for spin evolution, allow- 
ing for the formation of "asymmetric" BH binaries with 
01 7^ ^2- The model is not meant to rival the complex- 
ity of population-synthesis codes like StarTrack. Our 



goal is rather to isolate the physical ingredients that are 
specifically relevant to BH spin alignment. The model 
builds, when necessary (e.g. when computing the rem- 
nant masses resulting from gravitational collapse as a 
function of the progenitor masses, or in treating the CE 
phase) on results from StarTrack, and in Section IV we 



present a preliminary comparison of our conclusions with 
publicly available results from StarTrack. 



A. Length scales 

Before describing our astrophysical model, we review 
the length scales associated with the formation, inspiral, 
and merger of BH binaries. The well defined hierarchy in 
these length scales demonstrates the necessity of our joint 
analysis of astrophysics and PN evolution. GW emission 
[321 133| causes a binary with a semimajor axis less than 



Oh ^ 45 



q 



tew 



(1 + g)2 yiQio yrs / VlOM© 



(- 



M 



.1 1/4 



i?0 (4) 



to merge on a timescale tew less than the Hubble time 
— lO^^yrs. The astrophysical processes described in 
this Section, including mass transfer, SN explosion^ and 
CE evolution, are required to shrink the binary down to 
separations smaller than ajj. GW emission also circu- 
larizes the binary at separations comparable to ajj. PN 
spin-orbit couplings become important at much smaller 
separations 



OpNi 



10- 



GM 



10- 



M 



R. 



(5) 



below which they can lock binaries into resonant config- 
urations with well defined spin directions [B]. Previous 
studies of PN resonances for supermassive BHs j51-fTT| 
found that the effectiveness of resonance locking strongly 
depends on the orientation of the BH spins when the bi- 
nary reaches the separation apNi. The spin orientation is 
set by the binary's astrophysical formation history. Reso- 
nance locking can be important even at separations above 



flLIGO 



10" 



/ M 



V lOMe 



1/3 



GW 



V20Hz 



-2/3 



Re 



(6) 



at which the binary reaches the lower limit /gw — 
10 — 20 Hz of the Advanced LIGO/Virgo sensitivity band. 
The third-generation Einstein Telescope is expected to 
reach even lower frequencies of order /gw — IHz. Since 
these frequencies are well within the regime where PN 
resonances are important, a unified treatment of the as- 
trophysical initial conditions and of the subsequent PN 



^ Throughout the paper we will loosely use the term "supernova" 
to indicate the core collapse of massive stars, even when such 
events are not luminous. 
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evolution of the binary is essential to determining which 
spin configurations are most relevant for GW detectors. 
Such a treatment is the main goal of this work. 



B. Fiducial scenarios for binary evolution 

In this Section we describe how massive main-sequence 
binary stars evolve into BH binaries. Fig. [3] summarizes 
the critical stages of binary evolution in our model. To 
isolate the effects of spin orientation during the PN in- 
spiral of the BH binaries, we fix the final mass ratio to 
the typical value q — 0.8 [H]. To ensure that this final 
mass ratio is obtained, the initial stellar masses of the bi- 
naries must be fixed to (M^^,M^'J = {35Mq, 16.75Mq) 
in the SMR scenario, or {30Mq,24Mq) in the RMR sce- 
nario. Throughout the paper, we use a single prime to 
identify the initially more massive stellar progenitor or 
"primary" , and a double prime to denote the initially 
less massive progenitor or "secondary" . This choice of 
initial masses also fixes the total mass of our BH bina- 
ries to M = 13.5M0, quite close to the expected peak of 
the distribution for the total mass [12]. The mass of the 
stars is somewhat smaller than expected for the progen- 
itors of BHs of these masses because we have neglected 
stellar winds, that lead to considerable mass loss prior 
to BH formation. Table |l] provides numerical values for 
the masses and radii of both the primary and secondary 
throughout the evolution in both the SMR and RMR 
scenarios. Appendix |A 1| shows how this choice of initial 
masses leads to BHs of the desired final masses. 



SMR RMR 



SMR RMR 



35M0 
16.75M0 
3OM0 
8.5M0 
8Mq 
7.5M0 
6Mq 
17.97?0 
8128i?0 



3OM0 
24Mq 
35M0 
8M0 
8.5M0 
6M0 
7.5M0 
18.8i?0 
8787i?0 



^Sf 

R'c 

Rq 

R'a 

flnoCE 
fflmCE 



9.57Rq 

6.36i?0 
8.78Rq 
0.26Rq 
0.27 Rq 
36O8J?0 
3500^0 
6981J?0 
O.69J?0 



8.78i?0 
7.767?0 
9.577?0 
O.267?0 
O.277?0 
35OOi?0 
36O8i?0 
67587?0 
O.637?0 



TABLE I. Masses and length scales at various stages of the 
binary evolution in our SMR and RMR scenarios, as shown 
in Fig. [3] The only independent parameters are the main- 
sequence masses Mg; and Mg';, which have been tuned to 
study final BH binaries with mass ratio q — 0.8. The other 
values are defined in the main text, and they are obtained 
using the analytical prescriptions presented in Appendix [X] 



drawn from the distribution described in Appendix |A 2[ 
We assume that the spins of the primary S' and sec- 
ondary S" are initially alignecQ with the orbital angular 
momentum L. As the primary evolves, its envelope ex- 
pands until it fills its Roche lobe, initiating stable mass 
transfer to the secondary (phase b in Fig. [3|. The ef- 
ficiency of mass transfer is usually parametrized via a 
parameter fa G [0, I]: cf. Eq. (A9) of Appendix A3 We 



assume this mass transfer continues until the primary has 
depleted its hydrogen envelope, leaving behind a helium 
core of mass = 8.5Mq (M^ = 8Mq) in the SMR 
(RMR) scenario. Following [T^], we assume semiconser- 
vative mass transfer: the secondary accretes a fraction 
fa = 1/2 of the mass lost by the primary, growing to a 
mass M'^j = 30Mq (M^'^ = 35Mq) in the SMR (RMR) 
scenario at the end of the mass-transfer episode. In prin- 
ciple mass transfer should also change the orbital separa- 
tion, but we neglect this change as it is smaller than the 
width of the distribution of initial separations, as well 
as subsequent changes in the separation during the CE 
phase. 

Following the end of mass transfer, the primary ex- 
plodes in a SN (phase c in Fig. l3| producing a BH of 



M'bh 



7.5 Mq (M^ 



BH 



6 Mq ) in the SMR 



mass 

(RMR) scenario. For simplicity, in our simulations the 
spin of this newly born BH is assumed to be maximal 
{Xi = 1, * = 1 , 2) and aligned with its stellar progenitor. 
The SN ejecta are generally emitted asymmetrically, im- 
parting a recoil velocity to the BH which is generally a 
fraction of the typical recoil velocities for protoneutron 
stars: wbh ^ (l - /fb)i'pNS, wh ere /fb e [0, 1] is a "fall- 
back parameter" (cf. Appendix A 4). This recoil tilts the 



orbital plane by an angle 71, and changes the semimajor 
axis and eccentricity to ai and ei, respectively. These 
orbital changes depend on both the kick and the mass 
lost during the SN, as described in Appendix | A 5| 

After the SN explosion of the primary, the secondary 
evolves and expands. The primary raises tides on the 
swollen secondary, and dissipation may allow these tides 
to both circularize the orbit (so that the final eccentricity 
is et — 0) and align the spin S" of the secondary with 
the orbital angular momentum L, as shown in phase d of 
Fig. [Sj This tidal alignment is described in much greater 
detail in Appendix | A 6[ Given the uncertainty in the ef- 
ficiency of tidal alignment, we explore both extreme pos- 
sibilities: complete circularization and alignment of S" 
( "Tides" in Fig. [2| and no circularization and alignment 
at all ("No Tides" in Fig. 2]). As the secondary expands 
further, it fills its Roche lobe initiating a second phase 



The initial main-sequence stage of the evolution is 
shown as phase a in Fig. |3] Binaries are assumed to 
form on circulai|^ orbits with initial semimajor axes oq 



The initial eccentricity has minimal effect. In fact we have re- 



peated our calculations using an initially thermal distribution of 
eccentricities of the form /(e) = 2e, and we observed no signifi- 
cant difference in the final distribution of A'l' and 612 . 
The alignment of stellar spins in eclipsing binaries can be mea- 
sured through the Rossiter-McLaughlin effect 34 35 . Although 
many systems have aligned spins ^36 - 38 , there arc notable excep- 
tions [39j . We expect efficient tidal alignment in the progenitors 
of merging BH binaries, due to their small initial separations. 
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FIG. 3. A schematic representation of our model for BH binary formation and spin evolution. Empty circles represent stars, 
filled circles represent BHs. Phase (a) shows the initial main-sequence stellar binary. Mass transfer from the primary to the 
secondary (b) leads to a possible mass-ratio reversal. The first SN kick tilts the angle between the spins and the orbital plane 
(c). Tidal interactions can realign the stellar member of the binary (d). The second SN kick tilts the orbital plane again (e). 
Gravitational radiation shrinks and circularizes the binary before our explicit PN evolution begins (f). 



of mass transfer. However, unlike the first mass-transfer 
event, this second mass-transfer phase will be highly un- 
stable [lOHH] • Instead of being accreted by the primary, 
most of this gas will expand into a CE about both mem- 
bers of the binary. Energy will be transferred from the 
binary's orbit to the CE, ultimately unbinding it from 
the system. This energy loss shrinks the semimajor axis 
of the binary from ai to aicE, as shown in phase d of 
Fig. [Sj More details about CE evolution, including the 
relationship between ai and a ice, are provided in Ap- 
pendix [XTj After the secondary loses its hydrogen enve- 
lope, the remaining helium core has a mass = SMq 
(M^ = 8.5Mq) in the SMR (RMR) scenario, as hsted in 
Table HI 

After the end of CE evolution, the naked helium core 
of the secondary rapidly completes its stellar evolution 
and explodes as a SN, as shown in phase e of Fig. [3] 



This explosion produces a BH of mass Afgjj = GMq 
(Af^'jj = 7.5Mq) in the SMR (RMR) scenario, as listed 
in Table [T] We assume that this BH has a maximal spin 
that is aligned with the spin S" of its stellar progenitor, 
as we did for the primary. The SN leads to mass loss 
and a hydrodynamical recoil that change the semimajor 
axis and eccentricity of the binary to 02 and 62 , respec- 
tively. It also tilts the orbital plane by an angle 9 that 
can be calculated using the same procedure as given for 
the first SN in Appendix ] A 5| The tilt resulting from the 
second SN is generally much smaller than that from the 
first SN (8 ^ 71) due to the comparatively larger or- 
bital velocity following CE evolution. This tilt changes 
the angles between L and the spins S' and S" to 73 and 
72 , respectively. If tides efficiently align S" with L prior 
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to the second SN, these angles are given by 



quadrupole formula 



cos 72 — COS 7i cos 8 + cos (p' sin 71 sin 8 , (7) 
(tides) (8) 



cos 72 — COS 6 



where ip' is the angle between the projection of S' in 
the orbital plane before the SN and the projection of the 
change in L into this same initial orbital plane. If (p' is 
uniformly distributed (the direction of the SN kick of the 
secondary is uncorrelated with the spin of the primary), 
the second term on the right-hand side of Eq. ([T]) averages 
to zero, implying that 72 > 72 for most binariet]^ This 
is the mechanism for creating a binary BH population 
preferentially attracted to the A$ = ±180° family of 
spin-orbit resonances in the SMR scenario and the A$ — 
0° family of resonances in the RMR scenario, as shown 
in Fig. [I 

If tides are inefficient, 72 is instead given by 

cos 72 — cos ^ (no tides) (9) 
= cos 7i cos 8 — sin n7 sin 7i sin 8 , (10) 



where ^ is given by Eq. ( A21 1, and vj is the angle between 



the projection of S" into the orbital plane before the sec- 
ond SN and the separation vector between the members 
of the binary. If w is independent of ip' and uniformly 
distributee!^ the second term on the right-hand side of 
Eq. (10) also averages to zero, implying that 72 — 72'. 



The small scatter about this relation follows from the 
lesser influence of the second SN kick (8 <C 71), which 
implies that the identical first terms on the right-hand 
sides of Eqs. Q and ( 10 ) dominate over the differing 
second terms. This explains the lack of preference for 
either family of resonances in the "No Tides" scenario 
shown in Fig. [2j 

After the second SN, the BH binary is left in a non- 
relativistic orbit that gradually decays through the emis- 
sion of gravitational radiation, as shown in phase f of 
Fig. |3] We calculate how this orbital decay reduces 
the semimajor axis and eccentricity using the standard 



^ Well separated distributions of 72 and 72 require SN kick ve- 
locities that are comparable to the orbital velocity prior to the 
first SN, but much smaller than the orbital velocity before the 
second SN. Fortunately such kick velocities are well motivated, 
as described in Appendix I A 41 

® This assumption is well justified because the primary and sec- 
ondary spins process at different rates [Hi and VI2 given by 
Eqs. l |14| and | |15[ l below] and the precession timescale tp^e ~ 
is short compared to the time tg^ ~ 10® yrs between SN 
events. At lowest PN order, tpre ~ iLc('!'/c)~*, where tj^c = 
GM/c^ ~ 5 X 10~^{M/10Mq} s is the light-crossing time. At a 
separation a we have i)/c ~ 5 X 10~^ (M/IOMq^^^ (a/ R@)-'^^^ , 



dt _ 5 C^"3 

da ~ ~6ACPAP 



a3 {l + qY ,^ 2^7/2 A 73 2 



37 
96^ 



de 19 e, „ / 121 ,\ / 73 n 37 , 

24 96 



da 12a('~^"H'+ 304^' 



(11) 



(12) 



To an excellent approximation, the BH spins simply pre- 
cess about L during this stage of the evolution, leaving 73 
and 72 fixed to their values after the second SN. Once the 
semimajor axis reaches a value opNi = lOOOM (in units 
where G = c = 1), we integrate higher-order PN equa- 
tions of motion as described in Section |III| to carefully 
model how the orbit and spins evolve. We assume that 
radiation reaction circularizes the orbit (epN = 0) by the 
time we start integrating the higher-order PN equations 
describing the precessional dynamics of the BH binary. 
This assumption is fully justified, as we will show by ex- 
plicit integration in Section [III B| below. 



C. Synthetic black-hole binary populations 

In the previous Section, we presented fiducial scenarios 
for the formation of BH binaries characterized by three 
choices: 

i) stable mass transfer prior to the first SN can pre- 
serve (SMR) or reverse (RMR) the mass ratio of 
the binary; 

ii) hydrodynamic kicks generated by the SN can have a 
polar {Ob = 10°) or isotropic {9i, — 90°) distribution 
with respect to the exploding star's spin; 

iii) tides do or do not circularize the orbit and align the 
spin S" of the secondary with the orbital angular 
momentum L prior to the second SN. 

In this Section, we construct synthetic populations of 
BH binaries for the 8 different scenarios determined by 
the three binary choices listed above. To generate mem- 
bers of these synthetic populations, we perform Monte 
Carlo simulations^ in which random values determine 



i) the initial semimajor axis oq (Appendix A 2), 



ii) the magnitude and direction of the kick produced 
in the first SN (Appendix | A 4[), 



iii) the magnitude and direction of the kick produced 
in the second SN (Appendix A 4), 



so to 



■ 0.5 yr < tsN- 



^ We generated 10* binary progenitors to calculate the rates listed 
m Table which are therefore accurate to within ^ 0.01%. To 
avoid cluttering, we only show a subsample of 10* progenitors in 
the figures of this Section. 
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Kicks 


Tides 


Mass transfer 




/'(V \ 
l'mCE(%) 


!^SN2(%) 


m[/o) 


^'bh(7o) 


Isotropic 
Isotropic 


On 
On 


SMR 
RMR 


32.50 (80.50) 
32.55 (80.28) 


26.53 (12.24) 
34.86 (14.91) 


2.66 (0.51) 
2.97 (0.30) 


0.04 (0.00) 
0.04 (0.00) 


38.27 (6.74) 
29.59 (4.50) 


Isotropic 
Isotropic 


Off 
Un 


SMR 


32.50 (80.50) 

QO t^K /on OCA 

oz.oo (oU.Zo) 


26.53 (12.24) 
M.OO (i4.yij 


2.93 (0.60) 

o.Ui (^U.oOJ 


0.04 (0.01) 

U.U4 (^U.UUJ 


38.01 (6.65) 


Polar 
Polar 


On 
On 


SMR 
RMR 


31.84 (83.14) 
31.86 (82.97) 


26.68 (9.40) 
34.88 (12.10) 


3.29 (0.24) 
3.65 (0.24) 


0.01 (0.01) 

0.02 (0.00) 


38.18 (7.21) 
29.58 (4.70) 


Polar 
Polar 


Off 
Off 


SMR 
RMR 


31.81 (83.16) 
31.84 (82.98) 


26.65 (9.38) 
34.89 (12.09) 


3.35 (0.52) 
3.65 (0.33) 


0.03 (0.01) 
0.04 (0.00) 


38.15 (6.93) 
29.59 (4.60) 



TABLE II. Fraction of binaries v (in percentage) that satisfy the following conditions, each of which successively prevent the 
formation of a merging BH binary: i) are unbound by the first SN (i^sNi), ii) merge during the CE phase (t'mCE), iii) are 
unbound by the second SN {vs'N2), iv) do not merge within a Hubble time due to gravitational-radiation reaction (vb)- The 
final column is the fraction i^bh = 1 — iysNi + t'mCE + i^SN2 + vh) of all simulated binaries that form merging BH binaries. In 
parentheses we list the c orres ponding fractions if SN kicks are not suppressed by fallback (i.e. if we set /fb = rather than 
/fb = 0.8); see Appendix I A 41). 



iv) the angles tp' and vj specifying the directions of the 
spins S' and S" before the second SN (Section II B I, 



v) the angle A$ between the projections of the BH 
spins in the orbital plane at separation apNi- 

The angles vo, and A$ in items iv) and v) above 
are uniformly distributed in the range [0, 27r]. The syn- 
thetic populations generated in this procedure determine 
the initial conditions for the PN equations of motion de- 
scribed in Section Hill 

A binary-star system can fail to produce a merging BH 
binary for one of the foUoviring reasons: 

i) it is unbound by the first SN (ei > 1); 

ii) it merges during the CE evolution between the two 
SN {aicE < flmCE); 

iii) it is unbound by the second SN (e2 > 1); 

iv) the time t required for gravitational radiation to 
shrink the semimajor axis from a2 to apNi, found 
by solving the coupled PN equations (111 and ( 12 1, 
exceeds the Hubble time in — 10^° Gyr. 

Table [H] lists the fraction of simulated binaries vsNi, 
i^mCE, '^SN2, and vh that fail to produce merging BH 
binaries for reasons i) through iv) listed above, as well as 
the fraction j/bh = 1 - (^^SNi + t'mCE + '^sn2 + that 
do evolve into such binaries. 

The failure fractions indicate the relative importance 
of different physical phenomena. To emphasize the sen- 
sitivity of our results to the highly uncertain SN kicks, 
we also show how these fractions change when the BH 
kick wbh = (1 — /fb)wpNS fully equals that imparted to 
the protoneutron star (/fb — 0) rather than our canoni- 
cal choice (/fb = 0.8); see Appendix A 4 for more details. 



Stronger kicks unbind more binaries during the first SN, 
increasing t'sNi and thereby reducing the overall fraction 
vbh of binaries that survive to form BH binaries. This 
qualitatively agrees with results of detailed population- 
synthesis models; see models S, V8, and V9 in [33]. We 
adopt /fb = 0.8 in the remainder of the paper. 



Fig. |4] shows how the choices that define our fidu- 
cial scenarios affect whether SN kicks unbind the bina- 
ries. One result apparent from this plot (and supported 
by the failure fractions i^sN hsted in Table |ll]) is that 
the probability of unbinding the system depends only 
weakly on whether the SN kicks are isotropic or polar. 
This is consistent with the findings of fiS], which sug- 
gest mild sensitivity to 9b when the typical kick velocity 
''^BH ^ 50 km/s is small compared to the orbital veloc- 
ity vo ~ 2.4x lO3(M/3OM0)i/2(a/i?0)-i/2 km/s. Fig. [3] 
also shows the effect of tides on the fraction i/bh of BH bi- 
naries produced. In the absence of tidal dissipation ( "No 
Tides"), the binaries have nonzero eccentricity (e^ ^ 0) 
when the second SN occurs. Eq. (A17) shows that the fi- 
nal semimajor axis ay has additional dependence on the 
true anomaly tpi in this limit, broadening the distribu- 
tion of a/, as can be seen in the right panel of Fig. |4] 
The kicks can add coherently to the large orbital veloc- 
ities near pericenter of highly eccentric orbits, allowing 
binaries to become unbound even after CE evolution has 
reduced the semimajor axis: cf. the handful of light-gray 
(green) points with aicE ^ lOi?0 in the right panel of 
Fig. |4| . This increases the fraction ^'sn2 of binaries un- 
bound in the second SN when tides are "Off" in Table HIl 
The importance of CE evolution can be seen as well: vir- 
tually all binaries that fail to form a CE {aicE ^ lO'*i?0) 
are unbound by the second SN. Binaries bound tightly 
enough to survive the second SN almost always manage 
to merge through GW emission in less than a Hubble 
time (t-H < 1). 



III. POST-NEWTONIAN INSPIRAL 

A. Post-Newtonian equations of motion 

At large orbital separations, the dynamics of BH bi- 
naries in vacuum can be approximated by expanding the 
Einstein equations in a perturbative PN series, where the 
perturbative parameter is the ratio v/c of the orbital ve- 
locity to the speed of light. For historical reasons, one 
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FIG. 4. (Color online.) Scatter plot showing the change in the semimajor axis due to the first (left panel: ao — >■ ai) and second 
(right panel: aicE — >■ 12) SN. All plots refer to the SMR scenario, but the behavior in the RMR scenario is very similar. Darker 
(red) dots represent binaries that remain bound after each explosion, while lighter (green ) dot s correspond to binaries that are 
unbound. Dashed lines show the minimum post-SN semimajor axis a/,Min given by Eq. (A23I and the critical semimajor axis 
tSmCE given by Eq. ( A34 1 below which binaries merge during CE evolution. Kicks are too small to saturate the isotropic limit 
a/,Min for at < lO^i?©. 



usually says that a quantity is expanded up to fcPN or- 
der if all terms up to order (w/c)^'^ are retained. Follow- 
ing common practice in the general relativity literature, 
in this Section we will use geometrical units such that 
G = c= 1. 

The PN approximation can describe the evolution of 
stellar-mass binaries down to separations a ^ lOM (i.e. 
a - lO-''i?0 for a BH binary with M = IOMq), be- 
yond which fully nonlinear numerical simulations are 
needed [l^HiS] . GW detection templates depend on 
the binary parameters when the system enters the sen- 
sitivity band of the detectors, which is well into the 
regime where PN corrections are significant, but astro- 
physical models of BH evolution (as implemented e.g. 
in population-synthesis codes) have so far neglected all 
general-relativistic effects. The main goal of this Section 
is to show that solving the PN equations of motion is 
necessary to determine the orientation of BH spins when 
binaries enter the sensitivity band of GW detectors such 
as Advanced LIGO/ Virgo and the Einstein Telescope. 

The PN equations of motion and gravitational wave- 
forms for spinning BH binaries were derived by several 
authors (see e.g. [50 52 ). Our previous investigations 
of spin dynamics considered binaries on circular orbits; 



as shown in Section HIB below, this is an excellent ap- 
proximation for most binaries in our sample. They also 
included high-order PN terms such as the monopole- 
quadrupole interaction and the spin-spin self interactions 



[51-fTT|. that we report for completeness below. 

For circular orbits with radius a and orbital velocity 
V — (GAI/aY^^ , the "intrinsic" dynamics of a binary sys- 
tem depends on 10 variables: the two masses (toi, 7712), 
the spins Si and S2 and the direction of the orbital an- 
gular momentum L. At the PN order we consider both 
spin magnitudes and the mass ratio q remain fixed during 
the inspiral. This leaves 7 independent degrees of free- 
dom. Because BHs are vacuum solutions of the Einstein 
equations, there is only one physical scale in the problem 
(the total mass of the binary M). Rescaling all quanti- 
ties relative to the mass M, we are left with 6 "intrinsic" 
parameters. 

It is convenient to analyze the precessional dynamics in 
the frame where the direction of the orbital momentum 
L lies along the z-axis. If we take (say) the x-axis to be 
oriented along the projection of Si on the orbital plane 
(see Fig. 1 in [3]), we are effectively imposing 3 additional 
constraints just by our choice of the reference frame (2 
components of L and 1 component of Si are set equal to 
zero). Then the only 3 variables describing precessional 
dynamics are the angles 61, O2 and A$, as defined in 
Eqs. ([2| and The angle between the two spins 6*12 is 
related to the other independent variables as follows: 

cos 6*12 = sin 9i sin 6*2 cos A<I> -f cos 61 cos 92- (13) 

In summary, for any given binary with intrinsic param- 
eters (g, Xii X2), the precessional dynamics is encoded in 
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the variables {9i , 62 , A$) as functions of the orbital veloc- These variables can be evolved forward in time by inte- 
ity V or (equivalently) of the orbital frequency uj = /M. grating the following PN equations of motion: 



dt 



r2i X Si, 

= X S2, 



Mfl2 =r]v^ [2 



5/0 , 3 V 



2q 



2M2 



S2 - 3 (^L • S2 j L - 3g [L ■ Si j L 
Si-3(l-Si)l-^(l-S2H 



(14) 
(15) 



dL 



dt r]A'P dt 



(S1 + S2); 



(16) 



dv 32 71 „ 2 743 + 92477 o 

— = V < 1 — t; h t; 

5 M I 336 



V- /6 r^ / 113 mf 25 



34103 13661 59 2 , VX1X2 , -r^ua r \/(i t\ oA^ia a \ 
''1 + TZ'''! ^ 73 — ( '21(Si • L)(S2 • L) - 247(Si • S2) 



18144 2016 
2 



18 



48 



s.i:(=F)'(™(s.-t)-^^3 



i=l,2 

16447322263 16 



. 4159 + 158767; 



672 



1712 



7 



139708800 3 
4415 358675 



ilR + In 47;) + I TT" 

105 ^ ' ^ V 48 



451 



56198689 \ 541 



217728 )^ 896 



5605 
2592' 



91495 



77 H 77 

4032 6048 ' 1512 ' 



(17) 



where r/ = 17111712/ M'^ and je — 0.577 is Euler's constant. 
The leading terms in Eqs. (14)-(15), up to 0{v^) or 



2.5PN order, describe precessional motion about the di- 
rection of the orbital angular momentum L. We assumed 
that these terms dominated during the PN inspiral of 
the previous Section, allowing 72 and 72 to remain fixed 
at a > lOOOAf . Spin-orbit couplings appear at 3PN, 
and they are the reason for the existence of the resonant 
configurations [6 . From Eq. ( 16 1 we see that the direc- 



tion of the angular momentum evolves on a precessional 
timescale, while Eq. (171 implies that its magnitude de- 



creases on the (longerjradiation-reaction timescale due 
to GW emission. The leading (quadrupolar) order of 
Eq. (17 1 is equivalent to the circular limit of Eq. (11) 
when we recall that = M/a. 

Higher-order PN terms in the equations of motion were 
recently computed [53^ . We modified Eqs. ( 14 1-( 17 ) to in- 
clude these new terms, finding that they affect the late- 
time dynamics of individual binaries but have negligible 
influence on the statistical behavior of our samples. The 
robustness of these statistical properties under the inclu- 
sion of higher-order PN terms was already noted in |9l- 
[TT] . For completeness we retained the higher-order PN 
terms that will be reported in [33] in our Monte Carlo 
simulations, but we stress again that they have no ob- 



servable impact on our results. 

At a given separation a, Schnittman's resonant config- 
urations can be found by forcing the three vectors Si, 
S2 and L to lie in a plane (A$ = 0°,±180°) and by 
imposing the constraint that the second time derivative 
of cos ^12 vanish [B] . A one-parameter family of config- 
urations with A$ = 0° and 9i < 62 satisfies this reso- 
nant constraint, as does a second one-parameter family 
with A* = ±180° and 9i > 62- As a decreases due 
to GW emission, the curves determined by these one- 
parameter families change, sweeping through a large re- 
gion of the (01, 62) parameter space. The resonant con- 
straint evolves toward the diagonal ^i = ^2 as a — > 0. 
Individual resonant binaries move towards the diagonal 
in the {9i, ^2) plane along trajectories over which the 
projection So • L of the spin combination Sq defined in 
the effective-one-body model [51] . 



(l + q)Si + (l + q-l)S2, 



(18) 



is approximately constant (cf. Figs. 1 and 2 of j9j). Res- 
onant configurations with A$ = tend to align the two 
spins with each other, so that ^12 0° near merger. On 
the other hand, configurations with A$ = ±180° identi- 
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fied by their constant value of So • L evolve towards 



cosyi2 



(l + g)So-L 

(xi + qx2)M^ 



(19) 



B. Initial conditions for the PN evolution 

By construction, all of the merging BH binaries pro- 
duced in Section [n] have M — 13. SM©, q — 0.8, and 
Xi — X2 — ^- For this mass ratio and these spin mag- 
nitudes, binaries become attracted towards resonances 
("resonant locking") at separations a < lOOM Pre- 
vious studies suggest that the spin-orbit resonances re- 
main influential provided 9 ^ 0.4 and Xi ^ 0-5 [5HTT]. To 
be safe, we begin following binaries at an initial separa- 
tion apNi = lOOOM large enough so that we can neglect 
spin-spin coupling at greater separations [S]. Recall that 
the mass ratio was defined such that q = 1712/ mi < 1. 
In the SMR scenario, the primary yields the larger BH 
(Mgjj > Mgjj), so the angles are initialized to be 



h = I21 



n = 72 



(20) 



In the RMR case, the primary transfers so much mass 
to the secondary prior to the first SN that it actually 
produces the smaller BH (-A/gn < -^^bh)i implying that 
we must reverse our initialization: 



^2 =72 



(21) 



Although our decision to neglect spin-spin coupling for 
a > opNi allows us to initialize 9i in this manner, the 
lower-order spin-orbit coupling allows A<I> to evolve on 
the precessional timescale, which is short compared to 
the time it takes to inspiral from 02 to apNi- We can 
therefore choose A$ at opni to be uniformly distributed 
in the range [—180°, -1-180°]. Finally, since gravitational 
radiation is very efhcient at circu lari zing; the orbit [to 
leading order e cx a^^/^^; see Eq. (12)], we assume that 



all BH binaries have circularized by the time they reach 
apNi- We checked this assumption by numerically in- 
tegrating Eq. (12 1 from 02 to apNi after initializing it 



with the values 62 predicted following the second SN; the 
residual eccentricity at apNi was less than 10"'* for all 
BH binaries in our samples. 



Results 



We evolved 10^ BH binaries for each of the 8 different 



fiducial astrophysical scenarios described in Section II C 
from an initial separatiorj^ apNi = lOOOA/ to a final sep- 
aration apNf = lOM. This final separation roughly in- 
dicates where the PN approximation breaks down and 



* The a = lOOOAf snapshots in the figures of this Section are taken 
shortly after the beginning of the PN evolution. The angle A$ 



full numerical relativity becomes necessary [16H49]. To 
reduce the Poisson noise in the histograms of Fig. [2j we 
used larger samples of lO'* BH binaries. We integrated 
the PN equatio ns ([T4l )-(17l using a STEPPERDOPRS inte- 
grator in CH — h j55l , progressively refining the time steps 
at small separations (see [9] for further details). 

In Fig. [5j we show the evolution of the dynamical vari- 
ables (6*1, 6^2, A$) for both the SMR and RMR scenarios 
with efficient tides and isotropic kicks. As already an- 
ticipated in the introduction, efficient tidal interactions 
lead to spin orientations that are strongly affected by 
spin-orbit resonances. When binaries are brought close 
enough to resonant configurations by precessional mo- 
tion and gravitational-radiation reaction, they no longer 
precess freely through all values of A<I>, but instead os- 
cillate about the resonant configurations OH]. In the 
SMR scenario, the initial orientation of the spins is such 
that 9i > 02, and the binaries lock into resonances with 
A$ — ±180° [darker (red) points in Fig. [5]. In contrast, 
in the RMR scenario the initial spins have di < 62 and 
the binaries lock into resonances with A$ = 0° [lighter 
(green) points in Fig. [s] . Once the binaries are trapped 
near resonances, they evolve toward the diagonal in the 
(01,62) plane, as seen in the left panel of Fig. [sj This 
corresponds to 612 — > 0° for binaries near the A$ = 0° 
family of resonances (RMR scenario). As seen in the 
right panel of Fig. [sj there is a much broader range of fi- 
nal values for 612 in the SMR scenario, because these final 
values depend on the ini tial astrophysical distribution of 
So • L according to Eq. (19 1. 



Fig. [6] shows that spin-orbit resonances can have an 
even stronger effect on BH binaries when SN kicks are 
polar (aligned within 6*6 — 10° of the stellar spin [55]). 
As discussed in Appendix |A 5[ exactly polar kicks tilt the 
orbital plane by an angle Q given by Eq. ( |A24[ ), which 



can only attain a maximum value cos" -*(2/3f^ (where 
/3 = Mf /Mi is the ratio of the total binary mass before 
and after the SN) without unbinding the binary. For 
/3 ~ 0.9, as in our SMR and RMR scenarios, 6 < 40°, 
and kicks are rarely large enough even to saturate this 
limit. This explains the much narrower distribution of 



varies on the precessional timescale and can therefore change 
quite rapidly before the separation decreases appreciably on the 
longer inspiral timescale. The initial clustering in A$ visible 
in the top-right panels of Figs. [5] and [6] is not a resonant ef- 
fect, as the binaries continue to sweep through all values of 
at these large separations. It results instead from the different 
rates at which binaries in the SMR and RMR populations pre- 
cess, segregating the groups from each other during the first few 
precessional cycles. This behavior is better illustrated by the 
animations available online at the following URLs, which refer 
to efficient tides with isotropic kicks, efficient tides with polar 
kicks, inefficient tides with isotropic kicks, and inefficient tides 
with polar kicks, respectively: 

http: //www. phy . olemiss . edu/~berti/tides_isotr .gif 
http : / /www. phy . olemiss . edu/~berti/tides_polar . gif 
http : / /www. phy . olemiss . edu/~berti/notides_isotr . gif 
jhttp: //www. phy . olemiss . edu/-berti/notides_polar . gif 
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FIG. 5. (Color online.) Scatter plots of the PN inspiral of maximally spinning BH binaries with mass ratio q = 0.8 from 
an initial separation apNi just above lOOOM to a final separation apNf = lOM. The left panel shows this evolution in the 
(01,02) plane and the right panel shows the evolution in the (A$,^i2) plane. Darker (red) and lighter (green) dots refer to 
the SMR and RMR scenarios, respectively. The initial distribution for these Monte Carlo simulations was constructed from an 
astrophysical model with efficient tides and isotropic kicks. An animated version of this plot is available online at the URL: 
|http: //www.phy ■ oiemiss . edu/-berti/tides_isotr .gif] 



initial values of 9i in the left panel of Fig. |6] compared to 
Fig. [5j Binaries with these smaller initial misalignments 
are more easily captured into resonances, as can be seen 
from the near total segregation of the SMR and RMR 
populations in A<i> by the time the binaries reach apNf — 
lOM in the right panel of Fig. [6| 

In our model, two physical mechanisms are responsi- 
ble for changing BH spin orientations: SN kicks and tidal 
alignment. Both mechanisms are critical: kicks generate 
misalignments between the spins and the orbital angular 
momentum, but only tides can introduce the asymme- 
try between these misalignments that causes one family 
of spin-orbit resonances (the A$ = ±180° family in the 
SMR scenario, the A$ = 0° family in the RMR sce- 
nario) to be favored over the other. When tidal effects 
are removed, as shown in Figs. [7]and|8j BH binaries are 
formed with 9i ~ 62 on average. Being symmetric un- 
der exchange of the two BHs, the evolution in the SMR 
and RMR scenarios is almost identical. As expected, 
the binaries do not lock into resonant configurations, in- 
stead precessing freely during the whole inspiral. In the 
late stages of inspiral, the binaries tend to pile up at 
A$ = ±90°, i.e. they spend more time in configura- 
tions where the projections of the two spins on the or- 
bital plane are orthogonal to each other. Unlike the spin- 
orbit resonances, configurations with A"I> — ±90° are not 
steady-state solutions to the spin-evolution equations in 
the absence of radiation reaction [6' . The pile up at these 



configurations however is an essential complement to the 
spin-orbit resonances for preserving the well known result 
that initially isotropic spin distributions remain isotropic 
(see e.g. [13'). The physical origin of this phenomenon 
merits further investigation. 



IV. COMPARISON WITH POPULATION 
SYNTHESIS 

We have demonstrated that viable astrophysical for- 
mation channels can result in BH binaries that are 
strongly affected by spin-orbit resonances during the late 
PN portion of the inspiral but before the binary enters 
the GW detection band. Therefore PN resonances can af- 
fect the observed dynamics of precessing binaries. Even 
more interestingly, the distribution of the angles A$ and 
9i2 is a diagnostic tool to constrain some of the main 
physical mechanisms responsible for BH binary forma- 
tion (namely the efficiency of tides, and whether mass 
transfer can produce mass-ratio reversal). 

However, some caveats are in order. Even our limited 
exploration of the parameter space of BH binary forma- 
tion models has shown that the influence of PN reso- 
nances depends sensitively on highly uncertain factors, 
such as the magnitude and direction of SN kicks, or the 
mass ratio and semimajor axis of the binary at various 
stages of its evolution. In this Section, we argue that: (i) 
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FIG. 6. (Color online.) Scatter plots of the same quantities shown in Fig. [s] for an astrophysical model with efficient tides and 
polar kicks. For an animated version of this plot, see: |http: //www.phy . olemiss . edu/~berti/tides_polar .gif | 
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FIG. 7. Scatter plots of the same quantities shown in Fig.|5]for an astrophysical model with inefficient tides and isotropic kicks. 
For an animated version of this plot, see: |http : //www . phy . olemiss . edu/~bert i/not ides_isotr . gif | 



our fiducial scenarios are indeed representative of the pre- A. Is our fiducial scenario representative? 

dictions of more sophisticated population-synthesis mod- 
els (Section IV A); and (ii) as a consequence, observations 

of spin-orbit resonances through their GW signatures can In our study we chose to follow the evolution of 

provide valu a.ble in sight into BH binary formation chan- two binary progenitors in detail, using a specific for- 

nels (Section |IVB|. mation channel. The resulting BH binaries resemble at 

least qualitatively the low-mass BH binaries that can be 
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FIG. 8. Scatter plots of the same quantities shown in Fig. [5] for an astrophysical model with inefficient tides and polar kicks. 
For an animated version of this plot, see: h ttp : / /www . phy . olemiss . edu/~bert i/not ides_polar . gif | 



formed through a wide range of compact-object forma- 
tion scenarios at a range of metallicities: see e.g. [12]. 

An important assumption made in this study is that 
of negligible mass loss. Current calculations suggest that 
the progenitors of the most commonly detected BH bi- 
naries will in fact have low metallicity and strongly sup- 
pressed mass loss |12]. The advantage of our approach 
is that by neglecting mass loss and focusing on a pair 
of fiducial binaries we can perform a "controlled experi- 
ment" to highlight how different physical phenomena in- 
fluence the efficiency of PN resonance locking. Variations 
in the range of initial binary masses, wind mass loss and 
other mass transfer modes will affect the mass distribu- 
tion of the binaries and the initial distribution of the mis- 
alignment angles (^i ,^2): but not our main qualitative 
predictions, that should be rather robust. 

This study included what we believe to be the most 
important physical mechanisms that could trap binaries 
in resonant configurations, but it is certainly possible 
that additional ingredients overlooked in our model could 
complicate our simple interpretation of the results. For 
example, our argument relies on a universal and deter- 
ministic relationship between stellar masses and compact 
remnants. By contrast, some studies suggest that the 
relationship between the initial and final mass may de- 
pend sensitively on interior structure [57], rotation, or 
conceivably even stochastically on the specific turbulent 
realization just prior to explosion. As a concrete exam- 
ple, recent simulations of solar-metallicity SN explosions 
by Ugliano et al. [57] (including fallback) and O'Connor 
and Ott [5^ (neglecting fallback) have produced non- 
monotonic relationships between the progenitor and fi- 



nal BH masses. Likewise, our argument makes the sensi- 
ble assumption that BH spins are aligned with the spin 
of their stellar progenitor, but neutron star observations 
suggest that the protoneutron star's spin axis may be 
perturbed in a SN [55]. 

Our case studies of binary evolution omit by construc- 
tion many of the complexities present in more fully de- 
veloped population-synthesis models. The inclusion of 
additional physics presents interesting opportunities for 
a more detailed understanding of the connection between 
poorly constrained assumptions in population-synthesis 
models and GW observations. Some of the limitations 
we imposed on our model - and therefore, interesting op- 
portunities for follow-up studies - are listed below: (1) 
we follow the formation and evolution of only two pro- 
genitor binaries, rather than monitoring a distribution of 
masses; (2) we only consider maximally spinning BHs, 
while we should consider astrophysically motivated spin 
magnitude distributions; (3) we adopt very simple pre- 
scriptions for mass transfer and evolution, which have 
minimal feedback onto the structure and evolution of 
each star; (4) we employ an extreme "all or nothing" limit 
for tidal interactions; (5) we assume that BHs are kicked 
with a specific fraction of the overall SN kick strength; 
(6) we neglect stellar mass loss, magnetic braking and 
other phenomena that can occur in different formation 
scenarios. 

In summary: while our fiducial scenario provides a rep- 
resentative environment to explore the physics of PN res- 
onances, the specific mass distribution and the quantita- 
tive distribution of the misalignment angles at the be- 
ginning of the PN-driven inspiral will depend on detailed 
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binary-evolution physics which is neglected by construc- 
tion in our toy model. It will be interesting to initialize 
our Monte Carlo simulations using more comprehensive 
binary-evolution models that include a distribution of 
progenitor masses, track tidal backreaction on the spins 
and orbit, and model in more detail mass transfer and the 
modifications it introduces to core and stellar evolution. 



B. Observational payoff 

Let us provide a specific example to illustrate these un- 
certainties and their potential observational payoff. Our 
fiducial model assumed relatively low-mass BHs. These 
systems receive strong SN kicks (due to small fallback) 
and are more significantly influenced by CE contraction 
(because of the greater relative effect of the envelope 
binding energy). By contrast, more massive BHs in the 
Star Track sample will accrete a significantly higher frac- 
tion of their pre-SN mass, which drastically suppresses 
the typical kick magnitude. As a result, massive BH bi- 
naries can be expected to have BH spins more aligned 
with the orbital angular momentum. 

This sort of qualitative difference between low- and 
high-mass BH binaries presents an opportunity for GW 
detectors. The most easily measurable quantity in GW 
observations is the "chirp mass" Mchirp = rf^^M, where 
M = 1711+1712 is the total binary mass and 77 — mim2/M^ 
is the symmetric mass ratio (see e.g. [BDJ [5T]). There- 
fore, even though current simulations suggest that the 
detected sample will be dominated by high-mass, nearly 
aligned BH binaries, observations can clearly identify the 
low-mass sample, which should exhibit significant initial 
misalignment and more interesting precessional dynam- 
ics. Given the significant uncertainties in population- 
synthesis models, even upper limits on the spin-orbit mis- 
alignment for high-mass BH binaries would be extremely 
valuable, either to corroborate the expectation of strong 
alignment or to demonstrate the significance of SN kicks 
for high-mass BHs. 

Based on our prototype study, let us assume that each 
PN resonance is an unambiguous indicator of a specific 
fomation scenario: hypothetical GW measurements of 
angles A$ ~ ±180° mean efficient tides in the "standard 
mass ratio" (SMR) scenario; measurements of A<I> ~ 0° 
mean that mass reversal also occurred (RMR); finally, 
A$ ~ ±90° is an indication that tidal effects were ineffi- 
cient (cf. Fig. [1]). Under these assumptions, statistically 
significant measurements of A$ could directly identify 
how often each of the three formation channels (efficient 
tides, SMR; efficient tides, RMR; inefficient tides) occurs, 
for each binary mass. 

To illustrate how informative these measurements 
might be. Fig. [9] shows the relative number of merging 
binaries that undergo mass-ratio reversal as a function of 
chirp mass, as derived from the most recent StarTrack 
binary-evolution models |12j . The figure (which is meant 
to be purely illustrative) refers to Subvariation A of the 
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FIG. 9. Histograms of binaries that do (RMR) or do not 
(SMR) undergo mass-ratio reversal as a function of chirp 
mass according to the publicly available StarTrack data from 
http://www.syntheticuniverse.org/. For illustration, here 
we choose Subvariation A of the standard model, in the ter- 
minology of [12] . A comparison of the upper and lower panels 
shows the striking differences in the chirp-mass distribution 
resulting from different choices for the metallicity Z. 



"standard model" of Dominik et al. [T2]. Each panel 
shows the chirp-mass distribution of binaries that either 
do (RMR, dashed blue histograms) or do not (SMR, 
red solid histograms) undergo mass-ratio reversal. This 
distribution has characteristic "peaks" at specific values 
of the chirp mass at any given Z and it depends very 
strongly on composition, as we can see by comparing the 
two panels (which refer to Z/Zq = 1 and Z/Zq = 0.1, 
respectively). According to our model, measurements 
of A<i> for a large enough sample of binaries would al- 
low us to reconstruct the shape of these histograms as 
a function of chirp mass, potentially enabling new high- 
precision tests of binary evolution, above and beyond the 
information provided by the mass distribution alone. 

A preliminary assessment of the main features of 
population-synthesis models that could be probed by 
these measurements can be inferred from Table Hill There 
we list the overall fraction of BH binary systems that 
undergo mass-ratio reversal for several different binary- 
evolution scenarios explored in [12'. The most dramatic 
difference is due to composition: with few exceptions, 
models with solar composition {Z/Zq = 1) almost exclu- 
sively produce SMR binaries, while models with subso- 
lar composition {Z/Zq = 0.1) produce comparable pro- 
portions of SMR and RMR binaries. Furthermore there 
are clear trends in the ratio RMR/SMR as a function 
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Variation 


Subvariation A 


Subvariation B 


Subvariation A 


Subvariation B 








Z/Zq = 0.1 


Z/Zq = 0.1 


Z/Zq = 


1 


Z/Zq = 


1 








SMR 


RMR 


# 


SMR 


RMR 


# 


SMR 


RMR 


# 


SMR 


RMR 


# 







Standard 


63.2% 


36.8% 


32496 


66.8% 


33.2% 


17038 


91.9% 


8.1% 


10160 


92.9% 


7.1% 


8795 


1 




A = 0.01 


67.9% 


32.1% 


12368 


67.4% 


32.6% 


11401 


93.6% 


6.4% 


8171 


93.6% 


6.4% 


8171 


2 




A = 0.1 


62.7% 


37.3% 


27698 


65.2% 


34.8% 


16885 


88.9% 


11.1% 


11977 


92.1% 


7.9% 


8577 


3 




A = 1 


54.2% 


45.8% 


51806 


65.7% 


34.3% 


19415 


79.1% 


20.9% 


15820 


91.6% 


8.4% 


8442 


4 




A = 10 


50.1% 


49.9% 


50884 


62.9% 


37.1% 


17939 


73.2% 


26.8% 


14425 


91.6% 


8.4% 


8321 


5 




Mns = 3Mq 


62.5% 


37.5% 


32236 


66.2% 


33.8% 


16868 


91.6% 


8.4% 


9972 


92.8% 


7.2% 


8589 


6 




Mns = 2Mq 


62.3% 


37.7% 


32535 


65.9% 


34.1% 


16804 


91.5% 


8.5% 


9922 


92.5% 


7.5% 


8590 


7 




a = 132.5km/s 


58.2% 


41.8% 


36546 


63.1% 


36.9% 


18935 


88.9% 


11.1% 


11099 


89.6% 


10.4% 


9334 


8 




t'BH = UpNS (BHs) 


56.2% 


43.8% 


948 


72.5% 


27.5% 


207 


56.2% 


43.8% 


16 


0% 


100% 


2 


9 




DBH = (BHs) 


56.3% 


43.7% 


52832 


58.8% 


41.2% 


34569 


66.3% 


33.7% 


35267 


65.2% 


34.8% 


32547 


10: 


Delayed SN 


61.4% 


38.6% 


27310 


66.3% 


33.7% 


13841 


81.5% 


18.5% 


1032 


81.2% 


18.8% 


881 


11: 


Weak winds 


58.4% 


41.6% 


33872 


63.6% 


36.4% 


17765 


70.5% 


29.5% 


21786 


64.2% 


35.8% 


16182 



TABLE III. BH binary rates predicted by Star Track. RMR (SMR) is the percentage of binaries that do (not) experience 
mass-ratio reversal due to mass transfer; # indicates the total number of BH binaries in the sample. Each row refers to a 
different variation over the "standard model" . The variations illustrate the effect of changing one parameter (CE binding 
energy A, kick magnitude, etc.) with respect to the "best guesses" of the standard model. Each row also shows the effect of 
changing the metallicity Z and the Hertzsprung-gap donor prescription. In Subvariation A (B), binaries can (can not) survive 
a common-envelope event during the Hertzprung-gap phase; see [12] for details). 



of the envelope-binding-energy parameter A discussed in 
Appendix A 7 (compare variations 1 to 4) ; the strength 
of SN kicks (variations 8 and 9); and the amount of mass 
loss through winds (variation 11). These parameters are 
also well known to significantly influence the overall num- 
ber and mass distribution of merging binaries. 

In conclusion, while our model needs further testing 
and scrutiny against more complete population-synthesis 
calculations, it strongly indicates that GW measurements 
of A$ and 612 will provide a useful diagnostic of compact 
binary formation, complementary to the more familiar 
mass and spin measurements. In the next Section we 
conclude the paper with an overview of the challenges 
and rewards associated with these measurements. 



V. DISCUSSION 

Previous Monte Carlo studies of the spin-orbit reso- 
nances discovered by Schnittman [B] showed that spins 
tend to lock in a resonant plane if the binary has mass 
ratio q > 0.4 and the dimensionless spin magnitudes 
Xi ^ 0-5 as long as there is an initial asymmetry in the 
relative orientation of the spins with respect to the or- 
bital angular momentum, i.e. 61 7^ 62 [SHU]- 

In this work we built a toy model for BH binary for- 
mation focusing on the main physical ingredients that 
can produce such an asymmetry: SN kicks (that tilt the 
orbital plane every time a BH is formed), tidal inter- 
actions (that tend to realign the spin of the star that 
collapses later with the orbital angular momentum) and 
mass transfer (that can produce mass-ratio reversal, so 
that the heaviest BH corresponds to the lighter stellar 
progenitor). We showed that for stellar-mass compact 
objects formed at the endpoint of isolated binary evolu- 



tion the required conditions should ubiquitously occur. 

Perhaps more interestingly, we demonstrated that the 
angle A$ between the components of the BH spins in 
the plane orthogonal to the orbital angular momentum 
is in one-to-one correspondence with the BH formation 
channel that gave birth to the BH binary: if tides are 
efficient the PN evolution attracts the spins to the res- 
onant plane with A$ ~ 0° (A$ ~ ±180°) if mass re- 
versal does (does not) occur. When tidal effects are 
inefficient the spins precess freely, and they pile up at 
A$ — ±90° by the time the binary enters the band of 
advanced GW detectors. A preliminary comparison with 
detailed population-synthesis calculations suggests that 
the fraction of binaries in each family of resonant config- 
urations, both overall and as a function of (chirp) mass, 
should provide a highly informative diagnostic on some 
of the main uncertainties involved in binary-evolution 
physics (metallicity, binding energy of the CE, magni- 
tude of BH kicks). Measuring this fraction will require a 
large sample of BH mergers with sufficient signal-to-noise 
ratio, but hopefully such a sample will be obtainable by 
Advanced LIGO/Virgo after some years of operation at 
design sensitivity. 

Our initial study merits detailed follow-ups to assess (i) 
the potential accuracy of GW measurements of the pre- 
cessional parameters, and (ii) the information that can be 
extracted by comparison with population-synthesis mod- 
els. 

Detailed studies are required from the point of view of 
GW data analysis. We have assumed for simplicity that 
each PN resonance can be easily and unambigously dis- 
tinguished. In practice, accurate matched-filtering mea- 
surements of the angles A$ and 9i2 will need more work 
on the GW source- modeling front. Relevant issues here 
include the construction of gravitational-waveform tern- 
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plates adapted to resonant configurations, the develop- 
ment of specialized parameter-estimation strategies and 
the understanding of systematic (as opposed to statis- 
tical) errors for second- and third-generation detectors. 
Spin modulations are known to influence both the am- 
plitude and phase of the emitted radiation, and while 
there are several preliminary investigations of parameter 
estimation from spinning, precessing binaries, the direct 
measurement of parameters characterizing the spin-orbit 
resonances may require the inclusion of higher-order spin 
terms and/or higher harmonics in the waveform models. 

From an astrophysical standpoint, the observable dis- 
tribution of binary systems as they enter the detector 
band should be calculated (more realistically) by apply- 
ing our PN evolution to initial data derived from state-of- 
the-art binary population-synthesis models. In addition 
to corroborating our results, such a study will establish 
a comprehensive library of reference models that can be 
compared to observational data using Bayesian or other 
model-selection strategies: see e.g. [62H68] for previous 
efforts in this direction. Such a study is necessary also to 
make contact with other observables, such as as the rate 
and mass distribution of compact binaries. Only with a 
comprehensive and self-consistent set of predictions can 
we quantify how much the information provided by PN 
resonances complements information available through 
other observable quantities. 

In conclusion, the direct observation of resonant lock- 
ing will be challenging from a GW data-analysis stand- 
point. However the relatively transparent astrophysical 
interpretation of PN resonances makes such an investi- 
gation worthwhile. Even if only observationally accessi- 
ble for the loudest signals, these resonances will enable 
unique insights into the evolutionary channels that pro- 
duce merging compact binaries. In our opinion, more 
detailed studies of resonant locking in connection with 
population-synthesis models will offer a great observa- 
tional opportunity for GW astronomy. 
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Appendix A: Binary-evolution phenomenology 

Binary population synthesis relies on copious guidance 
from both observations and theory fSS|. Simulations of 
binary evolution that self-consistently account for stel- 
lar structure and mass transfer are computationally ex- 
pensive and depend on a wide variety of parameters 
[nni IZO]- Models that hope to generate astrophysically 
realistic binary populations must tabulate the results of 
these simulations and calibrate them against observations 
P71 [5U1 155] . Well developed algorithms exist to quickly 
generate large synthetic compact-binary populations sim- 
ilar to those produced in more expensive direct simula- 
tions |27l [3^ . In this Appendix, we use such population- 
synthesis models to justify and put into context the sim- 
ple procedure adopted in this paper. To further validate 
our model, we have also performed a handful of detailed 
binary-evolution calculations with the binary-stellar evo- 
lution BSE code by Hurley et al. [3U]. When adopting 
similar assumptions (i.e., low stellar mass-loss rates and 
large envelope binding energies), the BSE code produces 
qualitatively similar evolutionary scenarios to the proce- 
dure outlined in the text. The simple model and fiducial 
scenarios considered in this paper do not account for a 
thorough exploration of the parameter space, but they 
illustrate the essential physics and demonstrate that PN 
resonance locking can be the preferred outcome of astro- 
physically motivated BH binary formation channels. 



1. Single stellar evolution 

In this Section, we provide relevant information about 
the evolution of isolated stars. Main-sequence stars born 
with a mass Ms have a radius [7T] 



Rq 



1.33 



Ms 



0.555 



(Al) 



Massive, metal-rich main-sequence stars lose a substan- 
tial amount of mass via winds prior to going SN, but 
we neglect this mass loss for simplicity. The inclusion of 
wind mass loss in our model would reduce the mass of 
the hydrogen envelope available to be transferred to the 
secondary during the first mass-transfer event. While ne- 
glecting this mass loss quantitatively changes the binary 
evolution, we believe that it does not qualitatively alter 
our conclusions. Larger (and appropriately chosen) ini- 
tial stellar masses would lead to final BH binaries with 
masses comparable to those considered in our model even 
in the presence of winds. 

Stars with main-sequence masses in the range 25 Mq < 
Ms < 40 Mq evolve into supergiants with hclium-core 
masses well approximated by 



Mr 



O.lMs + 5Mg 



(A2) 
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(cf. top panel of Fig. 14 of '57') and radii [7^ 

(Mc/Mq)4-5 



Rg 

Rq 



4950 



1 + 4{Mc/Mq)^ 



0.5 . 



(A3) 



Once the hydrogen envelopes have been lost, the naked 
helium cores have radii 73 



log 



Rc 



-0.699 -I- 0.0557 log 



M, 



c 



Mr, 



0.172 



-2.5 



(A4) 

We neglect further evolution of the naked helium star 
before SN. For the large masses typical of BH progenitors, 
the naked helium cores have radiative envelopes and do 
not expand substantially during subsequent shell burning 
[?7l [7^ . After going SN, a main-sequence star leaves 
behind a BH of mass (bottom panel of Fig. 14 of [27 ) 



Mbh ^ O.SMs ~ 3Mc 



2. Initial semimajor axis 



(A5) 



The initial binary separation gq is drawn from a uni- 
form logarithmic distribution in the range [amin, Omax] 
[571 [75H77] . The upper limit Omax is chosen to ensure 
that the primary fills its Roche lobe during its super- 
giant phase, while the lower limit amin is chosen so that 
the secondary does not fill its Roche lobe after receiving 
mass from the primary. The Roche-lobe radius Rl of a 
star of mass TOq, in an orbit of semimajor axis a about a 
companion of mass to^ is [301 IZSj 



RL{a,ma,mp) 



0.49Q2/3 



0.6Q2/3 + in(i + Qi/3) 



(A6) 



where Q = ma/'mp, so the above limits are determined 
by the constraints 



RLiamiiK, M's„ M'si) = R[ 

i?L(amin,M^'^,M^) = ^" 



(A7) 
(A8) 



These limits are somewhat arbitrary, but different choices 
would not affect our main results. In fact, binaries that 
do not go through mass transfer (a > amax) are so widely 
separated that they are easily unbound by the first SN, 
while binaries where mass is transferred back to the pri- 
mary prior to this SN (a < a^in) will merge in the CE 
phase. These limits will therefore only affect the failure 
fractions presented in Table not the spin alignments 
of merging BH binaries. 



3. Stable mass transfer 

When a star fills its Roche lobe, gas will either be sta- 
bly transferred to its companion or form a CE about both 
members of the binary. Stable mass transfer is discussed 
in this Section of the Appendix, while CE evolution is 



discussed in Section A 7 In general, the stability of mass 
transfer depends on the donor star, the accreting star, 
and the mass ejected to infinity; as a first approxima- 
tion, stability criteria are usually implemented by simple 
thresholds on the binary mass ratio, as summarized in 
|42j and references therein. For our mass ratios, mass 
transfer from the primary to the secondary prior to the 
first SN will be stable, while mass transfer from the sec- 
ondary to the primary between the two SN events will 
lead to the formation of a CE. A fraction fa of the mass 
lost by the primary in the first mass-transfer event will 
be accreted by the secondary, increasing its mass to 



M'^s=M'<^,+fa{M's, 



M'c) 



(A9) 



Fully conservative mass transfer {fa — 1) preserves the 
total mass of the system, while all of the mass lost by the 
donor is ejected from the system in fully non-conservative 
mass transfer [fa — 0). We assume that stable mass 
transfer is semiconservative [fa = 1/2), in agreement 
with the standard model of Dominik et al. [T^] . Larger 
values of fa during this first mass-transfer event will tend 
to favor the RMR scenario over the SMR scenario. Since 
fa is directly tied to the fraction of binaries that undergo 
mass-ratio reversal in a given mass and mass-ratio range, 
our model suggests that it is potentially measurable via 
GW observations. For simplicity, we assume that tides 
and the mass transfer itself efficiently circularize the orbit 
(but see ^79, SO^ for recent investigations of mass transfer 
and circularization in eccentric binaries). 



4. Supernova kicks: magnitude and direction 

Following [7], we assume that asymmetric SN events 
impart hydrodynamical recoils to the newly formed pro- 
toneutron stars. We calibrate the magnitude of this pri- 
mordial kick using observed proper motions of young 
pulsars: each protoneutron star is kicked with a veloc- 
ity WpNs drawn from a single MaxwcUian with parameter 
a — 265 km/s [81 . A fraction /fb of this asymmetri- 
cally ejected material falls back onto the protoneutron 
star and is accreted as it collapses into a BH. This fall- 
back suppresses the magnitude of the final kick imparted 
to the BH to ubh — (1 — /fb)^'pNs; for BHs with masses 
-^BH = (6M0, 7.5M0), as in our fiducial scenarios, sim- 
ulations suggest /fb ~ 0.8 [SHIHSI. This BH kick distri- 
bution is consistent with the observed proper motions of 
galactic X-ray binaries hosting BHs (HH [85] . Although 
our results are not extremely sensitive to the precise mag- 
nitude of the BH kicks, the existence of such kicks is 
crucial to our model, as they are the only observation- 
ally well motivated mechanism to introduce misalignment 
between the compact binary spins and the orbital plane. 

We assume that the BH kicks are distributed in a dou- 
ble cone of opening angle 6b about the BH spin and con- 
sider two extreme scenarios: isotropic [Oi, = 90°) or polar 
[0}, = 10°) kicks. There is some observational [86 > 87 and 
theoretical [SSllHn] support for the polar model. However 
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vo X Li 




FIG. 10. Definitions of the angles used in Appendix |A 5| to 
study SN kicks. Before the SN, the members of the binary 
have a separation r and relative velocity vq. Mass loss and 
the SN kick Vk tilt the orbital angular momentum from Li to 
Lf while leaving the spin S unchanged. 



we examine both possibilities because this choice has a 
significant effect on the resulting binary orbits, as dis- 
cussed in Appendix |A 5| below. Our choice of 9b = 10° 
in the polar model was partly motivated by a compara- 
ble observed misalignment between the spin and proper 
motion of the Crab pulsar [S3] . 



5. Supernova kicks: influence on the orbit 

In this Section, we describe how SN kicks are imple- 
mented in our Monte Carlo calculations. The expressions 
provided below have been published previously either un- 
der more restrictive assumptions |25| or using different 
notation |30j . but we rederive them here for clarity and 
completeness. Each SN reduces the mass of the binary 
and imparts a kick to the newly produced compact rem- 
nant. We calculate how these effects change the Kep- 
lerian orbital elements by applying energy and angular- 
momentum conservation to the binary before and after 
the SN. As the duration of the SN explosion is short 
compared to the other stages of binary evolution, we 
assume that this orbital modification occurs instanta- 
neously. The definitions of the angles used in this Ap- 
pendix are illustrated in Fig. [TO] 

In the simulations reported in this paper we assume 



that the binary is on a circular orbit (ci = 0) and that the 
stellar spins are aligned with the orbital angular momen- 
tum — 0) when the first SN occurs (we have actually 
relaxed the circularity assumption in additional simula- 
tions not presented here, and we verified that this has 
a negligible impact on our conclusions). If tides are in- 
efficient, both of these simplifying assumptions will not 
hold, in general, for the second SN. Therefore here we 
present general expressions for the post-SN orbital ele- 
ments. These expression were first derived (to our knowl- 
edge) in [3D], but here we use notation similar to that of 
Kalogera [5S]. 

The binary separation r for a Keplerian orbit with 
initial semimajor axis and eccentricity can be ex- 
pressed as 



1 + Ci cos Ipi ' 



(AlO) 



where tpi is the true anomaly. Values for the true anomaly 
at the time of the SN are chosen by assuming that the 
explosion is equally likely to occur at any given time. The 
time t after the binary reaches pericenter is given by 



2tt , 



-t = E — Ci sin E, 



where 



P = 27r 



\GM, 



(All) 



(A12) 



is the period of a binary of total mass Mi. The eccentric 
anomaly E is related to the true anomaly ipi by 



cos V'i 



COSi? 



1 — e,; COS E 



(A13) 



We assume that t is uniformly distributed in the range 
[0, P] and derive the corresponding values of V'l from 
these relations. 

The direction of the kick velocity Vk is defined by a 
polar angle 9k and an azimuthal angle 4>k- Here 9k is 
the angle between Vj. and the pre-SN orbital velocity vq, 
and the axis defined by (pk = is chosen to be parallel 
to the orbital angular momentum L (see Fig. 10). The 



direction of the spin S of the collapsing star is specified 
by the angle 7^ between S and L and the angle w between 
the projection of S in the orbital plane and the separation 
f between the members of the binary. In terms of these 
angles, the angle 9j, between S and is given by 

cos 9p — — (sin 9k sin (j>k sin S + cos 9k cos S) cos ru sin 7^ 



+ (cos 9k sin S — sin 9k sin (f>k cos (5) sin m sin 7j 



sm tik cos cpk cos 7i 



(AM) 



where the angle S between the orbital velocity and line 
of separation is given in terms of the true anomaly by 



cos S 



Ci sin 



(1 + 2ei costpi + ef) 



1/2 



(A15) 
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In our Monte Carlo simulations, kick directions are drawn 
from uniform distributions in 0^, cos^fc, and w. Kicks 
confined to within an angle 6*^ of the stellar spin S are 
therefore implemented by repeated draws from this dis- 
tribution such that 



< 



Op>n 



(A16) 



The SN reduces the total mass of the binary from 
to Mf and changes the velocity of the exploding star from 
vq to vo + Vk- Applying energy and angular-momentum 



conservation to the binary before and after the SN, we 
find that the final semimajor axis a/ and eccentricity e/ 
are given by [5D] 



- 2Uk 



1 



1 + Ci cos ■i/'i 



2(/3 

1 + 2ej cos + e 



2\ 1/2 



l-ei 



cos 9k 



(A17) 



l-e 



1 + Uk 



1 



2(/3-l) 



- 2ej cos -f 

1 + Ci cos V'i 



1/2 



1 - e? 



cos Ok — 



2uk 



Bi sin tpi sin 9k sin 0^ 

1 + Ci cos t/jq 

1 + 2ei cos tpi + e ■ 



Ufe sm t^fc cos (pk 
l + Ci cos ^/'i ^ 



1/2 



COS 6*;^ 



(A18) 



where /3 = Mf/Mi and Wfc is the magnitude of the kick 
velocity normalized to the circular orbital velocity before 



the explosion, i.e. 



Uk = Vk 



GM,, 



(A19) 



If the right-hand side of Eq. (A18l is negative, Cf > 1 
and the SN has unbound the binary. For binaries that 
remain bound, the orbital plane is tilted by an angle Q 
such that 



cos 9 = 



1 + Uk 



1+Uk 



1/2 



1 + 2ei cos ipi + ef 
1 + 2ei cos -01 + ef 



1/2 



cos 9k — 



cos 9k 



Ci sin Tpi sin 9k sin 4 
I + Bi cos tpi 

Ci sin "0^ sin 9k sin 0^ 
\ + Ci cos V'i 



+ (1 - ef) 



Ufc sm cos 0fc 
1 + Ci cos ^Ai 



-1/2 



(A20) 



I 

and the angle between S and L is changed from 7^ to ^, where 
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cos^ 



1 + life 



1-ej 



Uk 



1 + 2ei cos + ef 



1/2 



COS 9k 



Ci sin ipi sin dk sin t/ffe 
1 + E j cos V'i 



cos7j 



1 + ei cos ip. 

1+Uk 



sm f/fc cos (pk sm 7^ sm vj 
1/2 



1 + 2ei cos ^Ai + e| 

2 ^ -1/2 



cos 9k 



Ci sin '0i sin 9k sin ( 
1 + ei cos ipi 



1 + 6; cos Ipi 



(A21) 



When S is aligned with L before the SN (7^ — 0), the 
tilt of the orbital plane equals the misalignment of the 
exploding star's spin = 0). 

The above expressions greatly simplify for initially cir- 
cular binaries. For example, the SN will disrupt the bi- 
nary if 

ul+2ukCOs9k + l-2l3 > (e, = 0) . (A22) 

The equations simplify even further if S and L are ini- 
tially aligned (7^ — 0), in which case exactly polar kicks 
are given by 9k ~ 7r/2, (f>k = 0. Exactly polar kicks larger 
than Uk > \J2ji — 1 always unbind the binary, while for 
isotropic kicks a bound tail of the distribution remains 
provided Mfc < 1 -I- \f2^. If kicks are confined to cones 
within an angle 9\j of L, the minimum final semimajor 
axis is 

a/,Min 



cos^ 



(e,; =0,7, = 0) 



(A23) 



2/3- 

exactly polar kicks {9\y — 0) can only increase the semi- 
major axis (fly Mill > fli): while isotropic kicks (0f, = 90°) 
can reduce the semimajor axis by at most a factor of 2 
(a/,Min = ai/2). 

Exactly polar kicks also add a significant component 
of angular momentum perpendicular to the initial orbital 
plane, leading to a strong spin tilt: 

cos e = , ^ , (e, = 0, 7, = 0) . (A24) 

However, the maximum tilt that polar kicks can produce 
while the binary remains bound is 

,-l^oaN-l/2 



(A25) 



By contrast, isotropic kicks can make the binary more 
tightly bound, allowing greater latitude for kicks to pro- 
duce bound systems with large spin misalignments. 

In the limit that the kick velocity is small compared to 
the orbital velocity (uk ^ 1), as should be the case for 
the second SN after CE evolution has reduced the binary 
separation, the tilt of the orbital plane is given by 



O = Uk sin 9k I cos < 



(e, =0,7, 



0). 

(A26) 



6. Tidal alignment 



As discussed in Sec. |IIB[ tidal dissipation can circular- 
ize the orbit of the binary and align the spin of the sec- 
ondary with the orbital angular momentum between the 
two SN explosions [Ml HOI HI] • A detailed treatment of 
the theory of tidal damping in massive stars is far beyond 
the scope of this paper, and relatively little data exists 
to calibrate these theoretical models if we wished to do 
so. We therefore only consider the two extreme possibili- 
ties: tides can either fully circularize the binary and align 
the spin of the secondary, or they are completely ineffi- 
cient. We provide order-of-magnitude estimates for tidal 
processes below; those interested in more details should 
consult one of the several excellent published reviews of 
tidal processes [69 l l9 H l92] . 

Tides should generally act on both members of the 
binary. However tidal effects on the BH can safely be 
ignored, given its small size. We therefore focus on tidal 
effects on the secondary between the two SN (phase d 
of the evolutionary scenario presented in Fig. [s]). If the 
secondary is fully convective, as expected for the core of 
a BH progenitor, convection causes internal damping on 



the viscous timescale t 



V ~ -/'\3MsRl/Ls)^/^ where 
Ms, Rs and Ls are the mass, radius and luminosity of 
the secondary, and 7 is a prefactor that depends on de- 
tails of the stellar structure [92]. The orbit evolves on 
the tidal-friction timescale 



9 (A/bh+M5)Mbh V^^; 



4 X 10"'^fc7 



g(i + 0)VioMQ; Vioi?, 



1 / Ms X'^V 



2/3 



VlOUe 



-1/3 / \ 8 



Rs 



yrs. 



(A27) 



where Mbh is the mass of the primary, Q = Mbu/Ms 
is the mass ratio at this stage of the evolution, and ^7 
is a constant of order unity depending on the internal 
structure of the star [HI]. Though the details depend on 
the initial stellar spin, tidal friction should synchronize 
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and align the spin of the secondary with the now circular 
orbit on this same short timescale |92) . 

The most notable feature of the tidal-friction timescale 
^tid given by Eq. (A27l is its extremely steep depen- 



dence on the ratio a/Rs- While the secondary remains 
on the main sequence with a radius given by Eq. (All, 



this ratio is typically 100 or greater for binaries that 
avoid merging during CE evolution. This implies that 
tidal alignment occurs on timescales much longer than 
the Hubble time Ih — 10^° yrs. However, once the sec- 
ondary evolves to fill its Roche lobe, its radius is given 
by Eq. (A6) and the ratio a/Rs becomes of order unity. 
This reduces the tidal-friction timescale well below typ- 
ical stellar-evolution timescales of a few million years 
(hydrogen-core burning) or even the briefer time 



tHG - 2.7 X 10^ 



/ Mc \' / Rc 



V IOM0 / V 10i?G 



/ is 



V lOUp 



yrs 
(A28) 



that the secondary spends on the Hertzsprung gap af- 
ter exhausting the hydrogen in its core (i.e., the Kelvin- 
Helmholtz timescale of the core). Since our fiducial sce- 
narios require the secondary to fill its Roche lobe prior to 
the second SN, one might expect tidal alignment to al- 
ways be efficient. Substantial uncertainties remain in the 
model however. Stars with partially radiative envelopes 
may have longer tidal- friction timescales [371 ES] i and the 
stellar core may not efficiently couple to its envelope, as 
suggested by recent Kepler observations of core-rotation 
rates Therefore, for completeness, we also explore 

the "extreme" alternative scenario of completely ineffi- 
cient tidal alignment. 

Being dissipative in nature, tidal interactions decrease 
the semimajor axis in addition to circularizing the or- 
bit. This change is small compared to that induced by 
CE evolution, as discussed in the next Section, and can 
therefore be neglected along with the orbital changes pro- 
duced by other phenomena (e.g. magnetic braking and 
mass transfer). 



7. Common-envelope evolution 

If the semimajor axis ai of the binary following the 
first SN is greater than anoCE, as determined from the 
constraint 



i?L(a„oCE,M^'/,M^H) 



R'n 



(A29) 



with Rq given by Eq. ( A3 1 , the secondary does not fill 
its Roche lobe and no CE evolution occurs. For smaller 
values of ai, we use conservation of energy to determine 
how much the binary's orbit shrinks during CE evolu- 
tion. The gravitational binding energy of the CE can be 
expressed as 



(A30) 



10= 
K 

\ 
o 




SMR 
RMR 



FIG. 11. The semimajor cixis aicE at the end of CE evolution 
as a function of its initial value ai in both the SMR and RMR 



scenarios. If ai > OnoCE, as given implicitly by Eq. ( A29 1, the 



secondary fails to fill its Roche lobe, no CE evolution occurs, 
and aicE = ai. If aicE < fflmCE, as given implicitly by 



Eq. (A34I, the helium core of the secondary fills its Roche 



lobe prior to the end of CE and the binary merges, failing 
to eventually form a BH binary. The nonlinear relationship 
between pre- and post-CE semimajor axes when ai < UnoCE 
results from the nontrivial dependence of the CE efficiency 
parameter A on ai, as given by Eq. |A31l. 



where Mgj is the mass of the secondary at the onset of 
CE evolution, Mgj — M'^ is the mass lost by the sec- 
ondary during this evolution, TZ = R^^ai, Algp M^^) 
is the Roche-lobe radius of the secondary at the onset 
of CE evolution, and A is a dimensionless parameter of 
order unity that depends on the mass and structure of 
the secondary, notably the location of the core-envelope 
boundary. Full stellar-evolution codes can be used to cal- 
culate the appropriate value of A for our BH progenitors 
[94H96] . We adopt an analytic fit to Fig. 3 of [l2], which 
summarizes the results of these calculations: 



A = ae-^^/«s+c 



(A31) 



where a = 0.358, b = 7.19 x lO'^, and c = 0.05. Conser- 
vation of energy during CE evolution implies 



2ai 

solving for aicE yields 



Eh = -- 



2a 



ICE 



O-ICE — fll 



Sf 



1 



2 ai Af^'f - 



A 7^ 



(A32) 



(A33) 
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If aicE is less than a^cE, as determined from the con- 
straint 



(A34) 



with i?p given by Eq. (A4), the hehum core of the sec- 
ondary itself fills its Roche lobe before the end of CE 
evolution. This leads to a prompt merger, preventing 
the eventual formation of a BH binary. Our final pre- 
scription for aicE as a function of Oi is shown in Fig. 
CE evolution is crucial to our model, shrinking the semi- 



major axis by a factor ~ 10^ and thereby allowing the 
eventual BH binary to merge in less than a Hubble time. 

Motivated by hydrodynamical simulations [HZl [35] and 
previous work on binary evolution, we neglect accretion 
onto the primary BH during CE evolution. These studies 
suggest that the BH accretes at substantially less than 
the Bondi-Hoyle rate during the evolution, accumulating 
< 0.1 Af0 in mass. Given this small change in mass, we 
are justified in ignoring any resulting changes in the BH 
spin [99] . As noted in Appendix A 1 we also neglect the 



expansion of naked helium stars, and therefore explicitly 
forbid a helium-star CE phase [73] . 
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